#include "quadrilateral.h"

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>

// Global variable... Bad bad Programmer, no donut for you!

boost::mt19937 mersenne_twister2;
boost::uniform_real<> real_dist2(0.0,1.0);
boost::variate_generator<boost::mt19937&, boost::uniform_real<> > random_number2(mersenne_twister2, real_dist2);

namespace t2c2 {
	
Quadrilateral::Quadrilateral(Material* _material,const Vector3D& _vertex1,const Vector3D& _vertex2,const Vector3D& _vertex3,const Vector3D& _vertex4) : Patch(_material) {
	vertices[0] = _vertex1;
	vertices[1] = _vertex2;
	vertices[2] = _vertex3;
	vertices[3] = _vertex4;
}

Quadrilateral::Quadrilateral(const Vector3D& _vertex1,const Vector3D& _vertex2,const Vector3D& _vertex3,const Vector3D& _vertex4) : Patch() {
	vertices[0] = _vertex1;
	vertices[1] = _vertex2;
	vertices[2] = _vertex3;
	vertices[3] = _vertex4;
}

Quadrilateral::Quadrilateral(void) : Patch() {

}

Quadrilateral::~Quadrilateral(void){

}

Vector3D Quadrilateral::get_vertex(unsigned short _index) const {
	return vertices[_index];
}

void Quadrilateral::set_vertex(unsigned short _index,const Vector3D& _vertex) {
	vertices[_index] = _vertex;
}

double Quadrilateral::area(void) const {
	const Vector3D edge1 = vertices[1] - vertices[0];
	const Vector3D edge2 = vertices[3] - vertices[0];
	const Vector3D edge3 = vertices[1] - vertices[2];
	const Vector3D edge4 = vertices[3] - vertices[2];

	return 0.5*(edge1.cross_product(edge2)).length() + 0.5*(edge3.cross_product(edge4)).length();
}

Vector3D Quadrilateral::get_normal(const Vector3D&) const {
	const Vector3D edge1 = vertices[1] - vertices[0];
	const Vector3D edge2 = vertices[3] - vertices[0];
	return (edge1.cross_product(edge2)).normalize();
}

Vector3D Quadrilateral::get_sample_point(void) const {
	double u;
	double v;
	u = random_number2();
	v = random_number2();

	//Barycentric coordinates
	Vector3D sample_point = vertices[0]*(1.0-u)*(1.0 -v) + vertices[1]*(1.0-u)*v + vertices[2]*u*v + vertices[3]*u*(1.0-v);
	
	return sample_point;
}

double Quadrilateral::intersection(const Ray& _ray) {
        //Ares Lagae and Philip Dutré ray-quadrilateral intersection

        std::numeric_limits<double> double_limit;
        const double infinity = double_limit.infinity();
        const double epsilon = double_limit.epsilon();
        //const double epsilon = 0.00001;

        Vector3D intersection_point;
        double u,v;

        Vector3D E_01 = vertices[1] - vertices[0];
        Vector3D E_03 = vertices[3] - vertices[0];
        Vector3D P = _ray.get_direction().cross_product(E_03);
        double det = E_01.dot_product(P);
        if (std::abs(det) < epsilon) return -1;
        double inv_det = 1.0/det;
        Vector3D T = _ray.get_origin() - vertices[0];
        double alpha = T.dot_product(P) * inv_det;
        if (alpha < 0.0) return -1;
        if (alpha > 1.0) return -1;
        Vector3D Q = T.cross_product(E_01);
        double beta = _ray.get_direction().dot_product(Q) * inv_det;
        if (beta < 0.0) return -1;
        if (beta > 1.0) return -1;

        if ((alpha + beta) > 1.0) {

    // Rejects rays that intersect the plane of Q either on the
    // left of the line V11V10 or on the right of the line V11V01.

                Vector3D E_23 = vertices[3] - vertices[2];
                Vector3D E_21 = vertices[1] - vertices[2];
                Vector3D P_prime = _ray.get_direction().cross_product(E_21);
                double det_prime = E_23.dot_product(P_prime);
                if (std::abs(det_prime) < epsilon) return -1;
                double inv_det_prime = 1.0 / det_prime;
                Vector3D T_prime = _ray.get_origin() - vertices[2];
                double alpha_prime = T_prime.dot_product(P_prime) * inv_det_prime;
                if (alpha_prime < 0.0) return -1;
                Vector3D Q_prime = T_prime.cross_product(E_23);
                double beta_prime = _ray.get_direction().dot_product(Q_prime) * inv_det_prime;
                if (beta_prime < 0.0) return -1;
        }
        // Compute the ray parameter of the intersection point, and
        // reject the ray if it does not hit Q.

        double t = E_03.dot_product(Q) * inv_det;
        if (t < 0.0) return -1;
        else { 
                intersection_point = _ray.get_origin() + _ray.get_direction()*t;
                return t; 
        }
        // Compute the barycentric coordinates of the fourth vertex.
        // These do not depend on the ray, and can be precomputed
        // and stored with the quadrilateral.  

        double alpha_11, beta_11;
        Vector3D E_02 = vertices[2] - vertices[0];
        Vector3D n = get_normal(Vector3D());

        if ((std::abs(n.x) >= std::abs(n.y)) && (std::abs(n.x) >= std::abs(n.z))) {
        alpha_11 = ((E_02.y * E_03.z) - (E_02.z * E_03.y)) / n.x;
        beta_11  = ((E_01.y * E_02.z) - (E_01.z * E_02.y)) / n.x;
        }
        else if ((std::abs(n.y) >= std::abs(n.x)) && (std::abs(n.y) >= std::abs(n.z))) {  
        alpha_11 = ((E_02.z * E_03.x) - (E_02.x * E_03.z)) / n.y;
        beta_11  = ((E_01.z * E_02.x) - (E_01.x * E_02.z)) / n.y;
        }
        else {
        alpha_11 = ((E_02.x * E_03.y) - (E_02.y * E_03.x)) / n.z;
        beta_11  = ((E_01.x * E_02.y) - (E_01.y * E_02.x)) / n.z;
        }
        if (std::abs(alpha_11 - 1.0) < epsilon) {    

        // Q is a trapezium.
                u = alpha;
                if (std::abs(beta_11 - 1.0) < epsilon) v = beta; // Q is a parallelogram.
                else v = beta / ((u * (beta_11 - 1.0)) + 1.0); // Q is a trapezium.
        }
        else if (std::abs(beta_11 - 1.0) < epsilon) {

        // Q is a trapezium.
                v = beta;
                u = alpha / ((v * (alpha_11 - 1.0)) + 1.0);
        }
        else {
                double A = 1.0 - beta_11;
                double B = (alpha * (beta_11 - 1.0))
                                - (beta * (alpha_11 - 1.0) - 1.0);
                double C = alpha;
                double D = (B * B) - (4.0 * A * C);
                double Q = -0.5 * (B + ((B < 0.0 ? -1.0 : 1.0)
                                * std::sqrt(D)));
                u = Q / A;
                if ((u < 0.0) || (u > 1.0)) u = C / Q;
                v = beta / ((u * (beta_11 - 1.0)) + 1.0); 
        }       
        intersection_point = _ray.get_origin() + _ray.get_direction()*t;
        return -1;

}



}
